#include <stdio.h>
#include <stdlib.h> 
#include <math.h>
#include <memory.h>
#include <omp.h>
#define TIMES 150
char *ABIN;
char *BBIN;
char *XBIN;
int X,Y,Z;
double THRESSHOLD;
int size;

void plus(double **a,double *b,double *r,int i,int j,int k)
{
	int YZ = Y*Z;
	int iYZ = i*Y*Z;
	int iYZ1 = (i-1+X)%X*Y*Z;
	int iYZ2 = (i+1)%X*Y*Z;
	int jZ = j*Z;
	int jZ1 = (j-1+Y)%Y*Z;
	int jZ2 = (j+1)%Y*Z;
	int k1 = (k-1+Z)%Z;
	int k2 = (k+1)%Z;
	if(j > 0 && j < Y-1)
		r[iYZ+jZ+k] = 	a[iYZ+jZ+k][0]*b[iYZ+jZ+k]+
						a[iYZ+jZ+k][1]*b[iYZ1+jZ+k]+
						a[iYZ+jZ+k][2]*b[iYZ2+jZ+k]+
						a[iYZ+jZ+k][3]*b[iYZ+jZ1+k]+
						a[iYZ+jZ+k][4]*b[iYZ+jZ2+k]+
						a[iYZ+jZ+k][5]*b[iYZ2+jZ2+k]+
						a[iYZ+jZ+k][6]*b[iYZ2+jZ1+k]+
						a[iYZ+jZ+k][7]*b[iYZ1+jZ1+k]+
						a[iYZ+jZ+k][8]*b[iYZ1+jZ2+k]+
						a[iYZ+jZ+k][9]*b[iYZ+jZ+k1]+
						a[iYZ+jZ+k][10]*b[iYZ1+jZ+k1]+
						a[iYZ+jZ+k][11]*b[iYZ2+jZ+k1]+
						a[iYZ+jZ+k][12]*b[iYZ+jZ1+k1]+
						a[iYZ+jZ+k][13]*b[iYZ+jZ2+k1]+
						a[iYZ+jZ+k][14]*b[iYZ+jZ+k2]+
						a[iYZ+jZ+k][15]*b[iYZ1+jZ+k2]+
						a[iYZ+jZ+k][16]*b[iYZ2+jZ+k2]+
						a[iYZ+jZ+k][17]*b[iYZ+jZ1+k2]+
						a[iYZ+jZ+k][18]*b[iYZ+jZ2+k2];
	else if(j == 0)
		r[iYZ+jZ+k] = 	a[iYZ+jZ+k][0]*b[iYZ+k]+
						a[iYZ+jZ+k][1]*b[iYZ1+k]+
						a[iYZ+jZ+k][2]*b[iYZ2+k]+
						a[iYZ+jZ+k][4]*b[iYZ+jZ2+k]+
						a[iYZ+jZ+k][5]*b[iYZ2+jZ2+k]+
						a[iYZ+jZ+k][8]*b[iYZ1+jZ2+k]+
						a[iYZ+jZ+k][9]*b[iYZ+jZ+k1]+														
						a[iYZ+jZ+k][10]*b[iYZ1+jZ+k1]+
						a[iYZ+jZ+k][11]*b[iYZ2+jZ+k1]+
						a[iYZ+jZ+k][13]*b[iYZ+jZ2+k1]+
						a[iYZ+jZ+k][14]*b[iYZ+jZ+k2]+
						a[iYZ+jZ+k][15]*b[iYZ1+jZ+k2]+
						a[iYZ+jZ+k][16]*b[iYZ2+jZ+k2]+
						a[iYZ+jZ+k][18]*b[iYZ+jZ2+k2]+						
						a[iYZ+jZ+k][3]*b[(i+X/2)%X*YZ+jZ+k]+
						a[iYZ+jZ+k][6]*b[(i+X/2-1)%X*YZ+jZ+k]+
						a[iYZ+jZ+k][7]*b[(i+X/2+1)%X*YZ+jZ+k]+
						a[iYZ+jZ+k][12]*b[(i+X/2)%X*YZ+jZ+k1]+
						a[iYZ+jZ+k][17]*b[(i+X/2)%X*YZ+jZ+k2];
	else if(j == Y-1)
		r[iYZ+jZ+k] = 	a[iYZ+jZ+k][0]*b[iYZ+jZ+k]+
						a[iYZ+jZ+k][1]*b[iYZ1+jZ+k]+
						a[iYZ+jZ+k][2]*b[iYZ2+jZ+k]+
						a[iYZ+jZ+k][3]*b[iYZ+jZ1+k]+
						a[iYZ+jZ+k][6]*b[iYZ2+jZ1+k]+
						a[iYZ+jZ+k][7]*b[iYZ1+jZ1+k]+
						a[iYZ+jZ+k][9]*b[iYZ+jZ+k1]+
						a[iYZ+jZ+k][10]*b[iYZ1+jZ+k1]+
						a[iYZ+jZ+k][11]*b[iYZ2+jZ+k1]+
						a[iYZ+jZ+k][12]*b[iYZ+jZ1+k1]+
						a[iYZ+jZ+k][14]*b[iYZ+jZ+k2]+
						a[iYZ+jZ+k][15]*b[iYZ1+jZ+k2]+
						a[iYZ+jZ+k][16]*b[iYZ2+jZ+k2]+
						a[iYZ+jZ+k][17]*b[iYZ+jZ1+k2]+
						a[iYZ+jZ+k][4]*b[(i+X/2)%X*YZ+jZ+k]+
						a[iYZ+jZ+k][5]*b[(i+X/2-1)%X*YZ+jZ+k]+
						a[iYZ+jZ+k][8]*b[(i+X/2+1)%X*YZ+jZ+k]+
						a[iYZ+jZ+k][13]*b[(i+X/2)%X*YZ+jZ+k1]+
						a[iYZ+jZ+k][18]*b[(i+X/2)%X*YZ+jZ+k2];
}

void amultiply(double **a,double *x,double *b)
{	
	int i;
	for(i = 0;i < X;i++)
	{
		int j,k;
		for(j = 0;j < Y;j++)
			for(k = 0;k < Z;k++)
				plus(a,x,b,i,j,k);
	}
}

int stop(double **a,double *x,double *b)
{	
	double t[X*Y*Z];
	amultiply(a,x,t);
	int i;
	double all = 0;
	
	for(i = 0;i < X*Y*Z;i++)
		all += pow(b[i]-t[i],2);
	
	all = sqrt(all);
	if(all < THRESSHOLD)
		return 1;
	return 0;
}

double dotproduct(double *x,double *y)
{
	int i;	
	double all = 0;
	for(i = 0;i < X*Y*Z;i++)
		all += x[i]*y[i];
	return all;
}

int main(int argc, char *argv[])
{	
	FILE *fr;
	int i;
	if(argc == 1)
	{
		ABIN = "case1/A.bin";
		BBIN = "case1/b.bin";
		XBIN = "case1/x0.bin";
		X = 360;
		Y = 180;
		Z = 38;
		THRESSHOLD = 1.878932301845534E-005;
	}
	else if(argc == 2)
	{
		int t;
		sscanf(argv[1],"%d",&t);
		if(t == 1)
		{
			ABIN = "case1/A.bin";
			BBIN = "case1/b.bin";
			XBIN = "case1/x0.bin";
			X = 360;
			Y = 180;
			Z = 38;
			THRESSHOLD = 1.878932301845534E-005;
		}
		else if(t == 2)
		{
			ABIN = "case2/A.bin";
			BBIN = "case2/b.bin";
			XBIN = "case2/x0.bin";
			X = 720;
			Y = 360;
			Z = 38;
			THRESSHOLD = 6.724779532126749E-005;
		}
		else
			exit(0);
	}
	else
		exit(0);
	
	size = omp_get_num_threads();
	double time1,time2,timeCacl = 0,timeInit = 0;
	time1 = omp_get_wtime();
	
	//Init
	double **a = malloc(sizeof(double*)*X*Y*Z);
	double *b = malloc(sizeof(double)*X*Y*Z);
	double *x = malloc(sizeof(double)*X*Y*Z);
	for(i = 0;i < X*Y*Z;i++)
	{
		a[i] = malloc(sizeof(double)*19);
		memset(a[i],0,sizeof(double)*19);
	}
	memset(b,0,sizeof(double)*X*Y*Z);
	memset(x,0,sizeof(double)*X*Y*Z);
	//IO
	fr = fopen(ABIN,"rb");
	for(i = 0;i < X*Y*Z;i++)
		if(fread(a[i],sizeof(double),19,fr) != 19)
			return 1;
	fclose(fr);
	fr = fopen(BBIN,"rb");
	for(i = 0;i < X*Y*Z;i++)
		if(fread(&b[i],sizeof(double),1,fr) != 1)
			return 1;
	fclose(fr);
	fr = fopen(XBIN,"rb");
	for(i = 0;i < X*Y*Z;i++)
		if(fread(&x[i],sizeof(double),1,fr) != 1)
			return 1;
	fclose(fr);
	
	//Init
	double *r = malloc(sizeof(double)*X*Y*Z);
	double *t = malloc(sizeof(double)*X*Y*Z);
	double *p = malloc(sizeof(double)*X*Y*Z);
	double *q = malloc(sizeof(double)*X*Y*Z);
	double alpha,beta,gama;
	amultiply(a,x,r);
	for(i = 0;i < X*Y*Z;i++)
		r[i] = b[i]-r[i];
	memcpy(p,r,sizeof(double)*X*Y*Z);
	amultiply(a,p,q);
	
	time2 = omp_get_wtime();
	timeInit = time2-time1;
	//Calc
	int round = 0;
	time1 = omp_get_wtime();
	while(round < TIMES)
	{
		gama = dotproduct(q,q);
		alpha = dotproduct(r,q)/gama;
		
		for(i = 0;i < X*Y*Z;i++)
			x[i] += alpha*p[i];
		
		if(stop(a,x,b))
			break;

		for(i = 0;i < X*Y*Z;i++)
			r[i] -= alpha*q[i];
		amultiply(a,r,t);
		beta = -dotproduct(t,q)/gama;
		
		for(i = 0;i < X*Y*Z;i++)
			p[i] = r[i]+beta*p[i];
		for(i = 0;i < X*Y*Z;i++)
			q[i] = t[i]+beta*q[i];
		round++;
	}
	time2 = omp_get_wtime();
	timeCacl = time2-time1;	
	printf("%d %lf %lf\n",round,timeInit,timeCacl);
	return 0;
}
